////////////////////////////////////////////////////////////////////////////////
//
// A class implementing a deme within a population, exchanging migrants with
// other demes
//
//
////////////////////////////////////////////////////////////////////////////////


#include "deme.h"
 
#include "cond_var.h"

#ifdef _SHOW_MUT_DISTR_
	#include <fstream.h>
	ofstream f_site_hits("mut_hits_distr.sum");
	static int count_simul=0;
#endif

//Debugging internal function
void SearchLine(TNode *anc, TNode *x, int *hits, int deme, int line, int depth);

//------------------------------------------------------------------------------

//Initialization of static members of TDeme
TMigrationMatrix* TDeme::Migrations=NULL;

//==============================================================================
//------------------------------------------------------------------------------
TDeme&
TDeme::operator=(const TDeme& D) {
   _growth_rate=D._growth_rate;
   _samples=D._samples;
   _deme_size=D._deme_size;
   _lineages=D._lineages;
   _lidum=D._lidum;
   _id=D._id;
   _min_deme_size=D._min_deme_size;
   _max_deme_size=D._max_deme_size;
   _proba_coal=D._proba_coal;
   Migrations=D.Migrations;
   _is_migration=D._is_migration;
   NodeList.de_allocate_list();
   int size=D.NodeList.ListSize;
   if (size) {                
      NodeList.allocate_list(size);
      for (int i=0; i<size; ++i)
         NodeList[i]=D.NodeList[i];
   }
   NodeList.de_allocate_list();
   size=D.OriginalNodes.ListSize;
   if (size) {
      OriginalNodes.allocate_list(size);
      for (int i=0; i<size; ++i)
         OriginalNodes[i]=D.OriginalNodes[i];
   }
   
   _polym_sites=D._polym_sites;
   return *this;
}

//------------------------------------------------------------------------------
int
TDeme::copy_node_pointers(const TTree& T, const int curr_lin, const long& max_lin) {
	 bool trace=0;
    int base=curr_lin;
    //This node list is dynamic and will change according to migrations and coalescences
    if(trace) cout<<"\nAllocating Nodelist("<<max_lin<<");";
    NodeList.allocate_list(max_lin);
    //Keeps a copy of the original nodes
    if(trace) cout<<"\nAllocating OriginalNodes("<<_samples.count()<<");";
    OriginalNodes.allocate_list(_samples.count()); //Christian: we needed more space!!!

    //Christian: Step through each sampling group in the deme
    for(int i=0; i<_samples.GetItemsInContainer(); i++) {
        for(int j=0; j<_samples[i].n(); j++, base++) {
           	//Stamp the node as originally belonging to the current deme
           	T(base)->set_deme(_id);
           	//Christian: Load the age and stat group too!
           	T(base)->set_time(_samples[i].gen());
           	T(base)->set_stat_group(_samples[i].stat_grp());
            //Then copy the node pointer to the local lists
           	NodeList[base-curr_lin]=T(base);
           	OriginalNodes[base-curr_lin]=T(base);
           	if(_samples[i].gen()>_oldest_samp)
               _oldest_samp=_samples[i].gen();
        }
    }
    //Fill up the node list with all empty pointers
    //so that no further size change is required during the generation
    //of the genealogy
    for (int i=_samples.count(); i<max_lin; ++i) {
        NodeList[i]=NULL;
    }
    if(trace) {
    	int i;
    	cout<<"\nNODE LIST\n------------";
    	for(i=0; i<max_lin; i++) cout<< *(NodeList[i]);

		cout<<"\n\nORIGINAL NODES\n----------------";
    	for(i=0; i<_samples.count(); i++) cout<< *(OriginalNodes[i]);
    }
    if(trace) printf("\nDone with copy_node_pointers");
    return 1;
};
//------------------------------------------------------------------------------
//Check and implements migrations
int          
TDeme::send_migrants(TDeme& sink, const int& i, const long& time) {
   bool trace=0; int n=0;
   if (!_lineages) return 0;
   //We will allow for multiple migrations per generation
   //Loro: !!! A faster routine could be made by first getting the number of migrants
   //from a binomial distribution
   for (int j=0; j<_lineages; j++) {
      //Christian: only extant lineages migrate, but if there is migration, allow what
      // used to be the last lineage a chance to migrate
      while (NodeList[j]->get_time()<=time && j<_lineages) 
        if (ran3(&_lidum)<=(*Migrations)(_id,i)) {
            migrate(sink,j);
            if(trace) printf("From %d to %d, Line %d (%d.%d) ",_id,sink._id,j, NodeList[j]->deme, NodeList[j]->node_number);
            n++;
        } else break;
   }
   if(trace && n>0) system("PAUSE");
   return 1;
};
//------------------------------------------------------------------------------
//Implements the coalescence process within population
//Christian: if 'hudson' is true then the next coalescent time may be appropriately chosen
//from an exponential distribution. The return value is the waiting time
//to the next coalescence. Otherwise, the return value is the number of lineages
//that coalesced in the current 'time'th generation
int TDeme::coalesce_lineages(const long& time, TTree& GeneTree, int pos, bool hudson) {
    bool trace=0;
    int extant_lineages=0, n_coals=0;
    if(trace) cout<<"\nGen" << time << ", Deme"<< _id << ":" << (int)_deme_size
        << "; " << _lineages << " lineages ";
    //Christian: count the number of lineages we're dealing with currently
    for(long i=0; i<_lineages; i++) {
        if(NodeList[(int)i]->time <= time)
            extant_lineages++;
    }
    if(trace) cout<<"(" << extant_lineages << " extant) ";
    //If there is only one EXTANT lineage left, do not even try to coalesce it...
    if (extant_lineages<2) return 0;
    
	if (!_deme_size) {
      cout 	<< "TDeme::coalesce_lineages() : deme size of population "
      		<< _id << " is zero!\n";
   	  return 0;
    }
    //prob = n(n-1)/(2N)
    double proba_coal;
    int skipgen=0;
    if(!hudson)
      proba_coal=_proba_coal*0.5*extant_lineages*(extant_lineages-1);
    else {
      proba_coal=1.0;
      skipgen=(int) (-2.0*log(ran3(&_lidum))/(_proba_coal*extant_lineages*(extant_lineages-1)));
    }
    
    if(trace) cout<< proba_coal;
    if (proba_coal==0.0) return 0;

   //Coalesce two lineages with probability proba
   while (ran3(&_lidum)<=proba_coal && extant_lineages>1) {
      //Pick two lineages at random for a coalescence
      //Christian: Don't allow samples older than 'time' to coalesce
      int pick1, pick2;
      if(trace) printf("\npicking [");
      long count=0;
      do {
         pick1=int(ran3(&_lidum)*_lineages);
         if(trace) printf("%d ",pick1);
         if(++count>1000) {printf("\n***BAD!***"); system("PAUSE"); return 0;}
      } while(NodeList[pick1]->time>time);
      if(trace) printf("] %d and [",pick1);
      count=0;
      do {
      	pick2=int(ran3(&_lidum)*_lineages);
         if(trace) printf("%d ",pick2);
         if(++count>1000) {printf("\n***BAD!***"); system("PAUSE"); return 0;}
      } while (pick2==pick1 || NodeList[pick2]->time>time);
      if(trace) printf("] %d",pick2);
      //Christian: end changes
      if(trace) cout<<"\n   *GOT ONE!!* Coalescing " << pick1 << " and " << pick2 << endl; 
      TNode *desc1=NodeList[pick1],
      		*desc2=NodeList[pick2],
            *ancestor=GeneTree(pos++);

      //Update node information
      ancestor->desc1=desc1;
      ancestor->desc2=desc2;
      ancestor->time=time+skipgen;
      desc1->ancestor=ancestor;
      desc2->ancestor=ancestor;

	  //Replace one coalescing node by the ancestor
      //and swap the second with the last node after updating the number
      //of remaining lineages 
      NodeList[pick1]=ancestor;
      NodeList[pick2]=NodeList[--_lineages];
      n_coals++;
      extant_lineages--;
      if(!hudson) proba_coal=_proba_coal*0.5*extant_lineages*(extant_lineages-1);
      else proba_coal=0.0;
   } if(trace) printf("\nNo more coals");
   if(hudson) {
      if(trace) printf("\n*Skipping %d generations",skipgen);
      return skipgen;
   }
   return n_coals;
};

//------------------------------------------------------------------------------
void
TDeme::reset() {
   //Restoring initial number of lineages
	_lineages=_samples.count();

    int size=NodeList.size();
    //Copying original nodes into dynamic list
	for (int i=0; i<_samples.count(); ++i) 
   	    NodeList[i]=OriginalNodes[i];
    //Filling up the remainint spaces in node list with empty pointers
    for (int i=_samples.count(); i<size; ++i)
        NodeList[i]=NULL;
    _min_deme_size=MIN_SIZE;
    _max_deme_size=MAX_SIZE;
    _proba_coal=1.0/(int)_deme_size;
}
//------------------------------------------------------------------------------
//A simple procedure to check whether there are migrations from this population
int
TDeme::check_migrations(){
   if (Migrations) {
      _is_migration=false;
   	int size=Migrations->size();
      for (int i=0; i<size; ++i)
      if (_id!=i) {
      	if ((*Migrations)(_id,i)>0.0) {
         	_is_migration=true;
            return 1;
         }
      }
      return 0;
   }
   return 0;
};
//------------------------------------------------------------------------------
void
TDeme::print_haplotypes_for_Arlequin(ostream & os, const int& n, const Mut_Type data_type) {

   char DNA_letters[4]={'A','G','C','T'};
   int len=OriginalNodes[0]->seq_length;

   os << "#Number of mutation hits per site\n"
   	<< "#Sites ";
   for (int j=0; j<len; ++j) {
      	os
            << setw(4)
            << setiosflags(ios::right)
         	<< j << " ";
      }      
   os << "\n";

   #ifdef _SHOW_MUT_DISTR_
   	f_site_hits << "Simulation #" << ++count_simul << "\t";
   #endif

   os << "#Hits  ";
   for (int j=0; j<len; ++j) {
      	os 
            << setw(4)
            << setiosflags(ios::right)
            << OriginalNodes[0]->hits[j] << " ";

         #ifdef _SHOW_MUT_DISTR_
      	f_site_hits
            << setiosflags(ios::right)
            << OriginalNodes[0]->hits[j] << "\t";
         #endif
      }
   os << "\n";

   #ifdef _SHOW_MUT_DISTR_
   f_site_hits  << "\n";
   #endif

   for (int i=0; i<_samples.count(); ++i) {
   	//Write id and frequency
   	os << n << "_" << i << " 1 ";
      TNode* CurNode=OriginalNodes[i];
      if (CurNode->sequence) {
      	int* cursite=CurNode->sequence->begin();
      	for (int j=0; j<len; ++j) {
         	if (data_type==DNA)
   					//os << DNA_letters[(*CurNode->sequence)[j]];
   					os << DNA_letters[*cursite++];
               else
   					//os << (*CurNode->sequence)[j];
                  //Loro_1_3_99 : In order to avoid negative numbers
               	if (data_type==MICROSAT) os << (10000+*cursite++) << " ";
               	else os << *cursite++;
         }
         os << "\n";
      }
   }
};
//------------------------------------------------------------------------------
void
TDeme::print_haplotypes_for_PAUP(ostream & os, const Mut_Type data_type) {
    bool trace=0;
    char DNA_letters[4]={'A','G','C','T'};
    int len=OriginalNodes[0]->seq_length;
    int n=_samples.count();
    if(trace) {
    	cout<<"\nPrinting "<< n << " haplotypes from deme "<<_id;
        for(int i=0;i<n;i++) cout<<"\n"<<i<<"\t"<<OriginalNodes[i];
    }
    for (int i=0; i<n; ++i) {
        TNode* CurNode=OriginalNodes[i];
        if(trace) cout<<*CurNode;
        if (CurNode->sequence) {    
            int* cursite=CurNode->sequence->begin();
			//Write id
            os << CurNode->node_number << "." << (CurNode->deme+1) << "\n";
            for (int j=0; j<len; ++j) {
     	      if (data_type==DNA) {
   			      //os << DNA_letters[(*CurNode->sequence)[j]];
   			      os << DNA_letters[*cursite++];
              } else if (data_type==MICROSAT) {
                    os << (10000+*cursite++) << " ";
                    //os << (*CurNode->sequence)[j];
                    //Loro_1_3_99 : In order to avoid negative numbers
       	      } else {
                os << *cursite++;
              }
            }
            os << "\n";
        }
    }
};
//------------------------------------------------------------------------------
//Christian 6/29/04 Oneliner, no breaks
ostream& operator<<(ostream& os, const TSamp& S) {
   char buf[80];
	sprintf(buf,"Size:%3d  Age:%6d  Deme:%2d  StatGrp:%2d",S._size,S._age,S._deme,S._stat_grp);
	os<<buf;
   return os;
}
//------------------------------------------------------------------------------
/*Christian 6/29/04: If the file contains ancient dates, then the file format
  is 'samps age deme'. This code is only appropriate for ancient sample files */
istream& operator >>(istream& is, TSamp& S) {
   BayesRead(&(S._size),&is,2);
   BayesRead(&(S._age),&is,2);
   BayesRead(&(S._deme),&is,1);
   BayesRead(&(S._stat_grp),&is,1);
   return is;
}

//------------------------------------------------------------------------------
ostream&
operator<<(ostream& os, const TDeme& D) {
//Christian: modified to include more than one sample per deme
   int n=D._samples.GetItemsInContainer();
	os << "\nDeme size   : " << (int)D._deme_size
   	<< "\nSamples     : " << n;
	for(int i=0;i<n;i++) os << "\n - " << D._samples[i];
   os << "\nGrowth rate : " << D._growth_rate
      << "\n";
    return os;
}

//------------------------------------------------------------------------------
ostream&
operator<<(ostream& os, const THistoricalEvent& HE) {
	os	  << "\nTime             : " << HE.time
   	  << "\nSource           : " << HE.source
        << "\nSink             : " << HE.sink
        << "\nMigrants         : " << HE.migrants
        << "\nNew size         : " << HE.new_deme_size
        << "\nNew growth rate  : " << HE.new_growth_rate
        << "\nNew migr. matrix : " << HE.MigMat
        << endl;
   return os;
}
//------------------------------------------------------------------------------
istream& operator>>(istream& is, THistoricalEvent& HE) {
	BayesRead(&HE.time,&is,2);
	BayesRead(&HE.source,&is,1);
   BayesRead(&HE.sink,&is,1);
   BayesRead(&HE.migrants,&is,3);
   BayesRead(&HE.new_deme_size,&is,3);
   BayesRead(&HE.new_growth_rate,&is,3);
   BayesRead(&HE.MigMat,&is,1);
   return is;
}
//------------------------------------------------------------------------------
//==============================================================================
//------------------------------------------------------------------------------
long
TDemeCollection::count_lineages() {
   long count=0;
   for (int i=0 ;i<_num_demes; ++i) {
   	count+=Demes[i]._lineages;
   }
   return count;
};
//------------------------------------------------------------------------------
 TDemeCollection&
 TDemeCollection::operator=(const TDemeCollection& DC) {
   _is_migration=DC._is_migration;
   _num_demes=DC._num_demes;
   _cur_mig_mat=DC._cur_mig_mat;
   if (Demes) delete[] Demes;
   if (DC.Demes) {
      Demes= new TDeme[DC._num_demes];
   	  for (int i=0; i<DC._num_demes; ++i) {
         Demes[i]=DC.Demes[i];;
      }
   }
   Events.Flush();
   if (DC.Events.GetItemsInContainer()) {
   	for (int i=0; i<DC.Events.GetItemsInContainer(); ++i) {
      	Events.Add(DC.Events[i]);
      }
   }
   MigrMatArray=DC.MigrMatArray;
   GeneTree=DC.GeneTree;
   if (CoalTimes) {
   	delete CoalTimes;
      CoalTimes=NULL;
   }
   if (DC.CoalTimes) {
      	try {
         	CoalTimes= new TMigrationMatrix();
        }
        catch (...) {
         	cout 	<< "TDeme::operator=(const TDeme& D): unable to allocate memory"
                  << endl;
            if (CoalTimes) {
      			delete CoalTimes;
      			CoalTimes=NULL;
      		}
            return *this;
        }
   	   *CoalTimes=*DC.CoalTimes;
   }
   if (PairDiff) {
   	delete PairDiff;
      PairDiff=NULL;
   }
   if (DC.PairDiff) {
      	try {
         	PairDiff= new TMigrationMatrix();
         }
         catch (...) {
         	cout 	<< "TDeme::operator=(const TDeme& D): unable to allocate memory"
                  << endl;
            if (PairDiff) {
      			delete PairDiff;
      			PairDiff=NULL;
      		}
            return *this;
         }
   	*PairDiff=*DC.PairDiff;
   }
   if (MinCoalTimes) {
   	delete MinCoalTimes;
      MinCoalTimes=NULL;
   }
   if (DC.MinCoalTimes) {
      	try {
         	MinCoalTimes= new TMigrationMatrix();
         }
         catch (...) {
         	cout 	<< "TDeme::operator=(const TDeme& D): unable to allocate memory"
                  << endl;
   			if (MinCoalTimes) {
      			delete MinCoalTimes;
      			MinCoalTimes=NULL;
      		}
            return *this;
         }
   	*MinCoalTimes=*DC.MinCoalTimes;
   }
   _polym_sites=DC._polym_sites;
   _tot_mut=DC._tot_mut;
   return *this;
}
//------------------------------------------------------------------------------
//Check for migrations and update coalescence probabilities within each population
int
TDemeCollection::check_migrations(){
   _is_migration=false;
   //Examine all populations in turns to know which
   //one experience migrations. !! Do not stop the checking after first migration found
   for (int i=0; i<_num_demes; ++i)  {
   	if (Demes[i].check_migrations()) {
      	_is_migration=true;
      }
   }
   if (_is_migration) return 1;
   return 0;
};
//------------------------------------------------------------------------------
int
TDemeCollection::check_growth(){
   for (int i=0; i<_num_demes; ++i)  {
   	if (Demes[i].growth()!=0.0) {
         return 1;
      }
   }
   return 0;
};
//------------------------------------------------------------------------------
int
TDemeCollection::create_demes(const TMigrMatArray* const MA) {
	MigrMatArray=(TMigrMatArray*)MA;   

	int size=(*MA)[0].size(); //Get number of demes from migration matrix
    TDeme aDeme;
    //Copies the pointer to the migration matrix in each deme
    aDeme.set_migration_matrix(&(*MigrMatArray)[0]);

    if (_num_demes!=size) {
   	    if (Demes) delete[] Demes;
        _num_demes=size;
        Demes= new TDeme[_num_demes];
    }

    for (int i=0; i<size; ++i) {
   	    //Set different ids for each deme
   	    aDeme.set_id(i);
   	    Demes[i]=aDeme;
    }
    return _num_demes;
};
//------------------------------------------------------------------------------
int
TDemeCollection::create_lineages(const TSampArray& SS) { 
//Christian 6/29/04: Modified for ancient sampling
   bool trace=0;
	initialize_samp_size(SS);

    //We need # of demes, not # of samples
    //int size=SS.GetItemsInContainer();
    //if (size!=_num_demes) return 0; 

    //Count total initial number of lineages (sum of samples, ancient and modern)
    long num_lineages=count_lineages();

    //Allocate room for the tree structure
    if (num_lineages) GeneTree.allocate_tree(num_lineages);
    else {
      cout << "TDemeCollection::create_lineages() : Tree is empty " << endl;
      return 0;
    }
   //Christian: GeneTree now contains all blank nodes. We need to load the
   // age information of the ancient samples. This is done in copy_node_pointers
   //Copy the pointers to the nodes of the genealogy in each deme
   for (int i=0, current_lineage=0; i<_num_demes; ++i) {
      Demes[i].copy_node_pointers(GeneTree,current_lineage,num_lineages);
      if(trace) printf("\n-Got the node pointers for Deme #%d",i);
      current_lineage+=Demes[i]._lineages;
   }
   if(trace) printf("\nDone creating lineages");
 	return 1;
};
//------------------------------------------------------------------------------
int
TDemeCollection::initialize_growth_rates(const TGrowthArray& GR) {
	int size=GR.GetItemsInContainer(); //Get number of demes from migration matrix
   if (size!=_num_demes) return 0;
  for (int i=0; i<size; ++i) {
   	Demes[i].set_growth_rate(GR[i]);
//      cout << i << "\t" << Demes[i].growth() << "\n";
   }
   return 1;
};
//------------------------------------------------------------------------------
int
TDemeCollection::initialize_deme_size(const TDemeSizeArray& DS) {
	int size=DS.GetItemsInContainer();
   if (size!=_num_demes) return 0;
   for (int i=0; i<size; ++i) {
   	Demes[i].set_deme_size(DS[i]);
   }
   return 1;
};
//------------------------------------------------------------------------------
//Loads the sample objects into the appropriate Deme
int
TDemeCollection::initialize_samp_size(const TSampArray& SS) { 
//Christian 6/29/04: modified to include ancient sample information
    bool trace=0;
    int size=SS.GetItemsInContainer();
    if(trace) cout<<"\n\nInitializing samples\n - # Samp: " << size << endl;
    //if (size!=_num_demes) return 0;
    for (int i=0; i<size; ++i) {
        //Demes[i].set_sample_size(SS[i]);
        if(trace) cout<<" - Sample " << i << ": Added to deme " << SS[i].deme();
        Demes[SS[i].deme()].AddSample(SS[i]);
        if(trace) cout<<" - Now " << Demes[SS[i].deme()];
        if(SS[i].stat_grp()>_num_stat_grps)
        	_num_stat_grps=SS[i].stat_grp();
    }
    return 1;
};
//------------------------------------------------------------------------------
int
TDemeCollection::initialize_events(const TEventArray& HE) {
	int size=HE.GetItemsInContainer();
   for (int i=0; i<size; ++i) {
   	Events.Add(HE[i]);
      //cout << Events[i] << endl;
   }
   return 1;
};
//------------------------------------------------------------------------------
int
TDemeCollection::implement_event(const int& cur_event) {
   THistoricalEvent& CurEvent=Events[cur_event];
   //Step 1: Send migrants from source to sink, fixed Christian 7/28/04
   // Check also that the sink population has a deme size >0
   if (CurEvent.migrants>0.0 && Demes[CurEvent.sink].deme_size()) {
      float x;
      long lidum=time(NULL);
      for (int i=0; i<Demes[CurEvent.source]._lineages; ++i) {
         x=ran3(&lidum);
         while(x<CurEvent.migrants && Demes[CurEvent.source].NodeList[i]->time<=CurEvent.time) {
            if(CurEvent.source!=CurEvent.sink) Demes[CurEvent.source].migrate(Demes[CurEvent.sink],i);
            else break;
            x=ran3(&lidum);
            if(i>=Demes[CurEvent.source]._lineages) break;
         }
      }
   }
   //Step 2: Resize sink deme size
   Demes[CurEvent.sink].linear_resize(CurEvent.new_deme_size);
   //Step 3: Readjust sink growth rate
   Demes[CurEvent.sink].set_growth_rate(CurEvent.new_growth_rate);
   //Step 4: Set new migration matrix; fixed Christian 7/28/04
   if (_cur_mig_mat!=CurEvent.MigMat) {
      //Copies the static pointer to the migration matrix in each deme
      for(int i=0; i<_num_demes; ++i) {
         Demes[i].set_migration_matrix(&(*MigrMatArray)[CurEvent.MigMat]);
         _cur_mig_mat=CurEvent.MigMat;
      }
   }
   return 1;
};
//------------------------------------------------------------------------------
int
TDemeCollection::build_tree() {
   bool trace=0;
   long num_lineages=count_lineages();
   long pos=2*GeneTree.sample_size()-num_lineages;
   long base=pos;
   int next_event=0, num_event=Events.GetItemsInContainer();
   int i, coal;
   //Loro_21_9_99 :
   //By default take the first migration matrix of the array: zero migration matrix
   // if no migration matrix is provided
   TDeme aDeme;
   //Copies the static pointer to the migration matrix in each deme
   aDeme.set_migration_matrix(&(*MigrMatArray)[0]);
   _cur_mig_mat=0;

   //Check that there are migrations in the genealogical process
   bool implement_migration=check_migrations(),
   		implement_growth=check_growth();
   		
   if(trace) {
      printf("\n\n**BUILDING TREE**\nInitial Conditions:\n");
      for(i=0; i<_num_demes; i++)
         cout<<Demes[i];
   }
   
   //We'll simulate a maximum of MAX_GENERATIONS generations
   for (long time=0, j, k; time<MAX_GENERATIONS && num_lineages>1; ++time) {
      if(trace) cout<<"\n\nGENERATION #"<<time<<" (" << num_lineages << ")\n-------------------";
      //if(time % 10000==0) cout<< "\n - Time="<<time <<" Lineages=" << num_lineages;
      //Step 1: Check for possible historical events
      if (next_event<num_event) {
        if(trace) cout<<"\n Checking historical events (next one: " << Events[next_event].time << ")";
      	while (Events[next_event].time==time) {
            if(trace) cout<<"\n - Doing one!";
         	implement_event(next_event);
         	//Check that the migration matrix has not changed
         	implement_migration=check_migrations();
         	//Check that the growth rates have not changed
         	implement_growth=check_growth();
            
            if(trace) {
               for(i=0; i<_num_demes; i++)
                  cout<<Demes[i];
               system("PAUSE");
            }
         	//Loro_21_9_99: To prevent the comparison of an unexisting event !  
            if (++next_event==num_event) break;
         } // end while event
      } //end if event

      //Step 2: Migration round
      if (implement_migration) {
        if(trace) cout<<"\n Migrating";
      	for (j=0; j<_num_demes; ++j) {
         	if (Demes[j].is_migration()) {
      			for (k=0; k<_num_demes; ++k) {
         			if (j!=k) Demes[j].send_migrants(Demes[k],k,time);
                }//end for
         	} //end if
      	} //end for
      } //end if

   	  //Step 3: Coalescence round; Christian: speed things up if the simulation
   	  //has reached a simple point by picking the next coalescent time from an
   	  //exponential distribution
      if(trace) cout<<"\n Coalescing";
      for (j=0; j<_num_demes; ++j) {
         if(next_event<num_event || implement_migration || implement_growth || Demes[j]._oldest_samp>time) {
         	coal=Demes[j].coalesce_lineages(time, GeneTree, pos, 0);
            if (coal) {
            	num_lineages-=coal;
               pos+=coal;
            }
         } else {
            if(Demes[j]._lineages>1) {
               time+=Demes[j].coalesce_lineages(time, GeneTree, pos, 1)-1; //returns gens skipped
               num_lineages--;
               pos++;
            }
         }
      } //end for Demes

   	//Step 4: Resize deme size after growth for next generation
      if (implement_growth) {
        if(trace) cout<<"\n Resizing demes";
        for (j=0; j<_num_demes; ++j) {
      	  if (Demes[j].growth()) {
            if (!Demes[j].exponential_resize()) {
	               cout << "TDemeCollection::build_tree(): size of deme "
               	     << j << "is zero !" << endl;
                   return 0;
            } //end if resize error
          } //end if Demes.growth()
          if(trace) cout<< "\n - Deme " << j << " Size:"  << Demes[j].deme_size();
        } //end for Demes
      } //end if implement_growth
      
      //Christian: Step 5: Verify that the original lineages are either still in the
      //deme collection, or descendents of nodes in the deme collection.
/* Take it on faith for now...
      TNode *x, *y;
      int n_fnd;
      for (j=0; j<_num_demes; ++j) {
         //if(trace) printf("\n\nOriginal Nodes from Deme %d\n",j);
         for(i=0;i<Demes[j].sample_size();i++) {
            x=Demes[j].OriginalNodes[i];
            //if(trace) printf("\n%d. %X",i,x);
            n_fnd=0;
            //search for it among the current lineages + their desc in all demes
            for(int k=0;k<_num_demes; ++k) {
               for(int l=0; l<Demes[k]._lineages; l++) {
                  y=Demes[k].NodeList[l];
                  SearchLine(y,x,&n_fnd, k, l, 0);
               }
            }
            if (n_fnd!=1) { printf("**ERROR** Node %d in deme %d has vanished!", i,j); system("PAUSE"); return 0;}
         }
      } //*/
    } //end for generations
    if (num_lineages>1) {
      cout << "TDemeCollection:: build_tree(): Coalescence process did not converge" << endl;
      return 0;
   }
   return 1;
};

void SearchLine(TNode *anc, TNode *x, int *hits, int deme, int line, int depth) {
   if (x==anc) {
      (*hits)++;
      /*if(depth) printf(" Gen%d desc of Line %d in Deme %d", depth, line, deme);
      else printf(" Line %d in Deme %d",line, deme);*/
   }
   if (anc->desc1) SearchLine(anc->desc1,x,hits,deme,line,depth+1);
   if (anc->desc2) SearchLine(anc->desc2,x,hits,deme,line,depth+1);
}


//------------------------------------------------------------------------------
void
TDemeCollection::reset(const TDemeSizeArray& SA, const TGrowthArray& GA) {
   initialize_deme_size(SA);
   initialize_growth_rates(GA);
   for (int i=0; i<_num_demes; ++i) {
   	Demes[i].reset();
   }
};
//------------------------------------------------------------------------------
//A procedure to sprinkle mutations once the tree branch lengths are known
// Modified declaration to add 12 parameters by Rumen Kostadinov
const long&
TDemeCollection::sprinkle_mutations(const float& mut_rate, const float& new_mut_rate, const long& time_of_mut_rate_change, 
				    const int& len,
				    const Mut_Type& mut_type, const float& gamma_par,
				    const float& mut_ratio, const float& prop_sites,
				    const float& trans_rate, 
				    const float& a_c, const float& a_g, const float& a_t,
				    const float& c_a, const float& c_g, const float& c_t,
				    const float& g_a, const float& g_c, const float& g_t,
				    const float& t_a, const float& t_c, const float& t_g,
				    const float& alpha, const float& beta,
				    const int& range_constr) {
  	//if(1) { GeneTree.count_descendants(); cout<<GeneTree; }

   //Get the root node
	TNode& root=GeneTree[2*GeneTree.sample_size()-2];

	//Set tree length to zero
   root.reset_tree_length();

   //Initialize hits per site to zero
   int size=root.hits.GetItemsInContainer(), to_add=len-size;
   //Loro_29_8_98 //Resize sequence length if needed
	for (int i=0; i<to_add; ++i){
         root.hits.Add(0);
   }
   for (int i=0; i<len; ++i)
   		root.hits[i]=0;
   long my_lidum=1L; //Any positive number will not reset the random number generator
   //Call a recursive procedure to add mutation according to a given mutation rate
   // Modified by Rumen Kostadinov
   // original: root.set_mut_rate(mut_rate);
   root.set_mut_rate(mut_rate, new_mut_rate, time_of_mut_rate_change);
   // Modified by Rumen Kostadinov
   // original:
   //_tot_mut=root.add_mutations(&my_lidum, 0, len, mut_type, gamma_par, trans_rate, range_constr);
   _tot_mut=root.add_mutations(&my_lidum, 0, len, mut_type, gamma_par, trans_rate, 
			       a_c, a_g, a_t,
			       c_a, c_g, c_t,
			       g_a, g_c, g_t,
			       t_a, t_c, t_g,
			       alpha, beta,
			       range_constr);
   return _tot_mut;
}
//------------------------------------------------------------------------------
extern ostream& operator<<(ostream& os,const TDemeCollection& DC) {
    int i=0, numdemes=DC.num_demes();
    os << "\nDEMES"
       << "\n======================";
	for (; i<numdemes; ++i) {
   	    os << "\nDeme #" << i+1
      	   << "\n----------\n"
           << DC.Demes[i]
           << "\n---------------------------------------\n";
   }
   int size=DC.MigrMatArray->GetItemsInContainer();
   for (int i=0; i<size; ++i) {
   	os << "\nMIGRATION MATRIX " << i
   		<< "\n=================\n"
      	<< (*DC.MigrMatArray)[i]
      	<< "\n";
   }
   os << "\nHISTORICAL EVENTS"
   	<< "\n=================\n";

   for (int i=0, numevents=DC.Events.GetItemsInContainer(); i<numevents; ++i) {
   	os << "\nEvent #" << (i+1)
      	<< "\n----------\n\n"
         << DC.Events[i]
         << "\n---------------------------------------\n";
	}
   return os;
};

//------------------------------------------------------------------------------
void
TDemeCollection::write_samples_to_Arlequin_file(ostream& os, const Mut_Type data_type) {
   for (int i=0; i<_num_demes; ++i)
   if (Demes[i].sample_size())  //Loro_20_9_99
   {
   	write_Arlequin_sample_header(i+1, Demes[i].sample_size(), os);
      Demes[i].print_haplotypes_for_Arlequin(os, i, data_type);
      write_Arlequin_sample_footer(os);
   }

}
//------------------------------------------------------------------------------
void
TDemeCollection::write_group_section_to_Arlequin_file(ostream& os) {
	os
   	  << "\n[[Structure]]\n"
      << "\n\tStructureName=\"Simulated data\""
      << "\n\tNbGroups=1"
      << "\n\tGroup={";
   for (int i=0; i<_num_demes; ++i) {
   	if (Demes[i].sample_size())   //Loro_20_9_99
   	os
      	<< "\n\t   \"Sample " << (i+1) << "\"";
   }
   os
   	<< "\n\t}" << endl;
};
//------------------------------------------------------------------------------
void
TDemeCollection::write_samples_to_PAUP_file(ostream& os, const Mut_Type data_type) {
   for (int i=0; i<_num_demes; ++i) {
		Demes[i].print_haplotypes_for_PAUP(os, data_type);
   }
}
//------------------------------------------------------------------------------
//Christian: The stat file!
void TDemeCollection::write_Statfile(ostream& os, const Mut_Type data_type) {
   bool trace=0;
	THapArray myhaps[_num_stat_grps+2];
	summary_stats gs;
	if(trace) cout<<"Writing stats";
	long ind, mrca=GeneTree[2*GeneTree.sample_size()-2].time;
	//Gather haptypes, make a combined group in myhaps[0]
	for(int curdeme=0; curdeme<_num_demes; curdeme++) {
		for(int curnode=0; curnode<Demes[curdeme].OriginalNodes.size(); curnode++) {
        	TNode *x=Demes[curdeme].OriginalNodes[curnode];
         T_HAP curhapc;
			if(trace) cout<<*x;
		   curhapc.hap=x->sequence;
		   //combined group
		   ind=myhaps[0].Find(curhapc);
			if(ind==INT_MAX) {
 			 	curhapc.n=1;
 			  	curhapc.is_private=1;
            curhapc.first_deme=curdeme;
		      myhaps[0].Add(curhapc);
            if(trace) cout<<"NEW"<<myhaps[0].GetItemsInContainer();
			} else {
				myhaps[0][ind].n++;
            if(trace) cout<<"MAT" << ind+1 << "." << myhaps[0][ind].n;
				if(myhaps[0][ind].first_deme!=curdeme) {
                myhaps[0][ind].is_private=0;
                if(trace) cout<<"*";
            }
			}
		   //individual stat group
		   T_HAP curhap;
		   curhap.hap=x->sequence;
			ind=myhaps[x->stat_group+1].Find(curhap);
			if(ind==INT_MAX) {
 			 	curhap.n=1;
 			  	curhap.is_private=1;
            curhap.first_deme=curdeme;
		      myhaps[x->stat_group+1].Add(curhap);
            if(trace) cout<<" NEW"<<myhaps[x->stat_group+1].GetItemsInContainer();
	      } else {
            myhaps[x->stat_group+1][ind].n++;
            if(trace) cout<<" MAT" << ind+1 << "." << myhaps[x->stat_group+1][ind].n;
				if(myhaps[x->stat_group+1][ind].first_deme!=curdeme) {
                myhaps[x->stat_group+1][ind].is_private=0;
                if(trace) cout<<"*";
            }
			}
			
		} //end for each node 
	} //end for each deme
	if(trace) { cout<< "\n--Finished loading haps--"; system("PAUSE"); }
   for(int curgroup=1; curgroup<=_num_stat_grps+1; curgroup++) {
      CalcStats(&(myhaps[curgroup]), &gs, data_type);
      if(data_type==DNA) {
         os << ","<<gs.num_haps << "," << gs.num_private_haps << "," << gs.seg_sites<<","
             << gs.pairws_dif<<","<< gs.hap_div <<","<<gs.nuc_div<<","<<gs.tajd <<",{";
          for(int i=0; i<gs.mismatch.GetItemsInContainer(); i++) os<<gs.mismatch[i]<<" ";
          os << "},";
      } else
         os << ","<<gs.num_haps << "," << gs.num_private_haps << "," << gs.hap_div<<","
             << gs.pairws_dif <<",";
      for(int curgroupy=curgroup+1; curgroupy<=_num_stat_grps+1; curgroupy++) {
         CalcInterStats(&(myhaps[curgroup]),&(myhaps[curgroupy]),&gs,data_type);
         if(data_type==DNA)
            os << ","<<gs.num_haps << "," << gs.num_private_haps << "," << gs.pairws_dif<<","
               << gs.hap_div <<","<<gs.nuc_div<<","<<gs.tajd <<",";
         else
            os << ","<<gs.num_haps << "," << gs.num_private_haps << ",";
      }
      myhaps[curgroup].Flush();
   }
	//combined stats
	if(trace) cout<<"\nCOMBO STATS\n";
	CalcStats(&(myhaps[0]),&gs, data_type);
	if(data_type==DNA) {
	   os << ","<<gs.num_haps << "," << gs.num_private_haps << "," << gs.seg_sites<<","
         << gs.pairws_dif<<","<< gs.hap_div <<","<<gs.nuc_div<<","<<gs.tajd <<",{";
      for(int i=0;i<gs.mismatch.GetItemsInContainer();i++) os<<gs.mismatch[i]<<" ";
      os << "},"<< mrca <<"\n";
   } else
  	   os << ","<<gs.num_haps << "," << gs.num_private_haps << "," << gs.hap_div<<","
         << gs.pairws_dif <<"\n"; 
	if(trace) { cout<<"\nDone.\n"; system("PAUSE"); }
}
//------------------------------------------------------------------------------

